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Abstract 

We consider the problem of computation and deformation of solutions of the com¬ 
plex Ginzburg-Landau equation (CGLE) with cubic nonlinearity in 1 -|-1 space-time 
dimension invariant under the action of the three-dimensional Lie group of symmetries 
A{x, t) —)■ A{x + (jA + t). From an initial set of such solutions, for a particular point 
in the parameter space of the CGLE, we obtain a sequence of new sets of invariant 
solutions via numerical continuation (homotopy) along paths in the parameter space. 
The solutions are obtained by solving an underdetermined system of nonlinear algebraic 
equations resulting from a spectral-Galerkin discretization of the CGLE in both space 
and time and the use of a path following method. 

The set of initial solutions led to distinct, new invariant solutions of the CGLE along 
the homotopy paths and in the final parameter region. The resulting solutions are un¬ 
stable, and have multiple modes and frequencies active in their spatial and temporal 
spectra, respectively. Symmetry gaining / breaking behavior associated mainly with 
the spatial reflection symmetry A{x, t) —>■ A{—x, t) of the CGLE was uncovered in the 
parameter regions traversed. 

Key Words: invariant solutions, relative periodic orbits, complex Ginzburg-Landau 
equation, continuous symmetries, numerical continuation 


1 Introduction 

We consider the problem of numerical computation and deformation of solutions of evolu¬ 
tionary partial differential equations (PDEs) fixed by the action of a Lie group F = Fi x M 
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of continuous symmetries of the PDEs, where M is the group of time translations and Fi 
is non-trivial. Within this context, such invariant solutions are also known as relative pe¬ 
riodic orbits or relative time-periodic solutions of an (autonomous) equivariant dynamical 
system. In this paper, we work with the complex Ginzburg-Landau equation with cubic 
nonlinearity in 1 + 1 space-time dimension, with Fi = (= x S^) - the two-torus. We 

note, however, that it should be straightforward to apply the methodology described in 
this paper to other evolutionary parameter-dependent PDEs invariant under the action of 
a group of continuous transformations. 

The complex Ginzburg-Landau equation (GGLE) is a widely studied PDE which has 
become a model problem for the study of nonlinear evolution equations exhibiting chaotic 
spatio-temporal dynamics, as well as being of interest in the context of pattern formation. 
It has applications in various fields, including fluid dynamics and superconductivity. (Eor 
details see, for example, [2, 17, 33, 38, 51] and references therein.) Eollowing [35], we 
consider here the GGLE with cubic nonlinearity in one spatial dimension, 

^ = 7271 + (1 + m)|^ - (1 + i^v)A\A\^ (1) 

with periodic boundary conditions 

A{x,t) = A{x + L^,t), (2) 

and spatial period Lx = 2 tt. The GGLE also appears in the literature in the form 
BA B'^ A 

-^ = A +{l + iv)-^ - {l + ipL)A\A\^, A{x,t) = A{x + L,t), (3) 

but note that with a change of variables 

Lx (Lx \ . Lx 

x — )■ — X, ^ A — A 

1j \ 1j J 1j 

one obtains equation (1), with R = (L/Lj,)^. Thus we adopt the formulation (1) without 
loss of generality and, henceforth, when we refer to the CGLE we mean equation (1) with 
the boundary conditions (2) unless otherwise noted. 

Equation (1) describes the time evolution of a complex-valued field A{x,t)- The pa¬ 
rameters R, ly, and /r in the equation are real. When R > 0 there is, in general, nontrivial 
spatio-temporal behavior and this is therefore the region of interest. The parameters u 
and /i are measures of the linear and nonlinear dispersion, respectively [2, 33]. It is known 
that in the form (1) the CGLE generates a continuous semiflow on a variety of spaces 
[5, 16, 33, 49]. 
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As will be discussed in detail in Section 2, the CGLE has a three-parameter group 
G = X M of continuous symmetries generated by space-time translations and a rotation 
of the complex field A. Thus, we focus our study on invariant solutions of the CGLE that 
satisfy 

A(x, t) = e“^A(x + S,t + T) 

for some (cp, S) G and T > 0. The interest here is on invariant solutions of the CGLE 
having multiple frequencies active in their temporal spectrum, not on single-frequency so¬ 
lutions A{x,t) = B{x)e^‘^^ (see, for example, [15, 24, 26]) or generalized traveling waves 
(also called coherent structures) A{x,t) = p{x — where p and cj) are real¬ 

valued functions and oj is some frequency (see, for example, [2, 9, 38, 50]). The class 
of single-frequency and generalized traveling wave solutions have been studied more ex¬ 
tensively than the multiple-frequency class, wherein lies our focus. The desired invariant 
solutions are sought from those of a system of nonlinear algebraic equations resulting from 
the use of a spectral-Galerkin discretization of the CGLE in both space and time. Hence, 
the method used does not fall in the (more commonly used) class of multiple-shooting al¬ 
gorithms for computing periodic orbits. Details concerning the discretization are provided 
in Section 3.1. We note that the CGLE is also invariant under the action of the discrete 
group of transformations A{x,t) A{—x,t) and thus solutions of the CGLE may also be 
hxed by this Z 2 -symmetry. While it is not uncommon in numerical continuation studies to 
center on solutions fixed by the Z 2 -symmetry (for example, even solutions), we make no 
such restriction here in order to work with a more general set of solutions. 

One should view the space of solutions and the space of G-orbits, or moduli space of 
solutions of the CGLE, as being fibered over the space of parameters The aim 

of the present study is to acquire a more global view of (a part of) the (fibered) space 
of G-orbits of the CGLE and its structure as we move around the point {R, i', p) in the 
CGLE parameter space. To this end, we start with a set of previously computed invariant 
solutions of the CGLE for a particular point (i?, v, p) in the parameter space of the CGLE 
[35] , and employ numerical continuation to carry solutions of this initial set (all of which are 
unstable) into solutions in a regime with different (fixed) parameter values. In particular, 
the new parameter region of interest was sought by increasing the value of R. However, 
at times we had to venture into a subspace (of the CGLE parameter space) with varying 
V and p as well. Taking the presence of positive Lyapunov exponents for typical (that 
is, non-invariant) solutions as an indication of chaotic dynamics [44], both the initial and 
final parameter regions in our study exhibit chaotic behavior. Specifically, non-invariant 
solutions in the initial and final parameter region have, respectively, 5 and 16 positive 
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Lyapunov exponents.^ This provides another motivation for conducting this study, which 
is to evaluate the potential benefits of using numerical continuation (on problems with 
a large number of unknowns) to continue multiple (distinct) unstable invariant solutions 
from one chaotic regime into another (with the aim, again, of ending with multiple (distinct) 
unstable invariant solutions in the final region), which is a topic of current interest [10]. 
(More generally, the problem of computing multiple unstable invariant solutions in a fixed 
parameter region is of interest, and the development/evaluation of different methods for 
this purpose is a subject of current research. See [35, 12, 10] and references therein for 
further details.) 

In more detail, we use a path following, or homotopy, method [43] to solve the un¬ 
derdetermined system of nonlinear algebraic equations resulting from the discretization of 
the CGLE, with the initial guesses for the nonlinear equations solver provided by the set 
of unstable solutions from [35]. Since the initial guesses are solutions of the PDE, it is 
then natural (and typical) to consider the parameters of the PDE as continuation (homo¬ 
topy) parameters. Clearly each successful step in the numerical continuation will yield an 
invariant solution of the PDE in the parameter regions traversed, and bifurcations along 
each continuation path may be encountered. In the present study, bifurcating behavior 
associated with symmetry breaking or gaining was frequently uncovered, which revealed a 
complex and interesting structure of the moduli space of G-orbits. We describe in detail in 
Sections 2 and 4 the additional symmetries that are being gained or broken at particular 
values of the parameters of the CGLE. In other words, we identified a set of special points in 
the CGLE parameter space at which the structure of the space of G-orbits changes abruptly 
(with a gain or loss of symmetry). This amounts to a kind of “phase transition” in the 
moduli space, and is one interesting aspect to consider as a focus for a further detailed 
investigation. 

The CGLE parameter values considered in [35] were chosen because for them the system 
exhibits temporal chaos, yet the value of R was such that a relatively small number of spatial 
modes could capture the dynamics of the system. The resulting problem sizes made the use 
of solvers for systems of nonlinear algebraic equations with a wider region of convergence 
than Newton’s method [13, 43], for example, the Levenberg-Marquardt method [32, 37, 41], 
amenable. In the absence of very close initial guesses, use of the aforementioned type of 
methods resulted in a favorable success rate, while use of Newton’s method was not a viable 
option [35]. Increasing the value of the parameter i?, though, leads to an increase in the 

^Lyapunov exponents for typical (that is, non-invariant) solutions were computed using the technique 
from [~, 8]. 
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number of unknowns required to capture well the dynamics of the system, which in turn 
eventually makes use of the aforementioned successful methods impractical or unfeasible 
due to memory and other computing requirements, particularly on personal workstations. 
However, with the availability of very close initial guesses, Newton’s method, which is 
typically used in conjunction with path following methods, could still be a suitable option. 

To the best of our knowledge, the way in which the Newton step is computed here 
(described in Section 3.2.1) is new. Although the particular approach is conceptually simple, 
that is where its value lies: it allowed us to compute an accurate Newton step quickly 
and efficiently and made the solution of a computationally challenging problem with a 
large number of unknowns (up to 32,260 unknowns in the present study) practical without 
the need of a cluster or supercomputer. We expand on these and other aspects of the 
numerical methodology employed, including the use of freely available packages to aid in 
our implementation of the required software, in Section 3. 

A detailed account of results from the current study is provided later on in Section 4. 
We note here that the set of solutions used as starting points in the path following method 
correspond to the first 15 solutions listed in the Appendix from [35]; these were selected 
simply to follow the order listed in said Appendix. The (unstable) solutions in this initial 
set were carried from the point {R, u, /i) = (16, —7, 5) in the CGLE parameter space to the 
regime with parameter values = (100,-7,5). The number of unknowns to solve 

for ranged between 4,000 and 32,260. Both the number of 15 solutions from [35] and the 
final point {R,u,^) = (100,-7,5) in the CGLE parameter space were chosen because we 
deemed them to be sufficient to help us gain insight into the structure of the solution space 
as well as to allow us to evaluate the potential for success of the proposed approach for 
computing multiple invariant solutions in fixed parameter regions of a dynamical system 
which exhibits chaotic behavior. 

The initial set of solutions led to distinct, new invariant solutions of the CGLE along 
the continuation paths and in the final parameter region. The resulting solutions are 
unstable, and have multiple modes and frequencies active in their spatial and temporal 
spectra, respectively. The fact that the computed solutions are unstable suggests that they 
may belong to the set of (infinitely many) unstable periodic orbits embedded in chaotic 
attractors [12, 10, 35]. This direction, by itself, is certainly very interesting to pursue in a 
future study of the dynamics of the CGLE. 

Previous studies on the CGLE which employ numerical continuation include [48], where 
bifurcations from a stable rotating wave to two-tori (of the aforementioned generalized 
traveling wave class) were identified. Values of L,,. = 1 and R < 180 were considered for 
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the CGLE (l)-(2) in the study [48], giving L ~ 13.42 for the maximum length of the 
spatial period in the formulation (3) of the CGLE. In comparison, the values = 27r and 
R < 100 in our study yield a maximum value of L ~ 62.83 in the formulation (3). The 
values of u and /r used in [48] are different from those used in the current study, but in 
both cases the parameters belong to the Benjamin-Eeir unstable region 1 + < 0 [51]. 

A different study [9] considers the aforementioned traveling waves solutions of the form 
A(x,t) = p{x — where p and cj) are real-valued functions and ui is some 

single frequency, in which case the CGLE reduces to a system of three coupled ordinary 
differential equations (ODEs) in z = x — vt. Continuation (on v and p,) was performed on 
the system of three ODEs for different values of L up to 512 and different chaotic regions 
were classified. 

Other studies can be found in [40], where transition to chaos from a limit cycle of the 
CGLE is investigated, [27], in which the bifurcation structure and asymptotic dynamics of 
even solutions of the CGLE are analyzed, and [34], where numerical continuation is used 
to study the dynamics of the CGLE in heteroclinic cycles, with a focus on invariant Z 2 - 
subspaces. A numerical study on solutions fixed by the spatial reflection, or Z 2 , symmetry 
of the CGLE and the stability or instability of such solutions with respect to perturbations 
that break said symmetry appears in [3], where parameter values of R = 1.05,16,36 are 
considered, and a spatial period of = 27r was used (the latter being the same as in 
the current study). The subsequent study [4] considers symmetry-breaking perturbations 
for solutions hxed by the spatial translation symmetry of the CGLE, for several (discrete) 
values of the parameter R in the range [4.2, 80]. 

Among more recent work, although not concerning the CGLE, we note [10], which 
deals with the numerical computation of invariant solutions embedded in a turbulent two- 
dimensional flow, with the aim of using periodic orbit theory to reconstruct turbulence 
statistics. Numerical continuation is considered as one of the alternatives for computing 
multiple invariant solutions of the two-dimensional Navier-Stokes equations for several fixed 
values of the Reynolds number. Among other aspects, the study in [10] documents the 
challenges arising when attempting to numerically compute multiple invariant solutions in 
hxed parameter regions (displaying chaotic or turbulent behavior), particularly when the 
problem involves large numbers of unknowns (potentially on the order of 10^ to 10®). 

The reader may consult [1, 14, 29, 46, 54] and references therein for a general back¬ 
ground on numerical continuation methods and software. In addition, dehnitions and 
results concerning symmetry groups and differential equations with symmetries can be 
found in [18, 22, 39, 42], and properties of relative equilibria and relative periodic orbits in 
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[19, 30, 31, 52, 53]. Finally, although outside the scope of this paper, the interested reader 
is referred to [12, 11, 55, 10] and references therein for treatments on the use of periodic 
orbits in the context of the study of chaotic dynamical systems, in which the computation 
of multiple unstable invariant solutions in fixed parameter regions is of much relevance. We 
proceed now with details concerning the problem addressed in the present study. 

2 Invariant Solutions of the CGLE and their Properties 

The CGLE has a number of well known symmetries that are central to its behavior [2]. In 
particular, equations (l)-(2) have a three-parameter group 

G = X M (4) 

of continuous symmetries generated by space-time translations x^x + a,t^t + T and 
a rotation A —)• A of the complex field A{x,t), in addition to being invariant under the 
action of the discrete group of transformations A{x, t) —)• A(—x, t) of spatial reflections. In 
other words, if A{x,t) is a solution of equations (I)-(2), then so are 


A{x, t), 

(5) 

A{x + a, t), 

(6) 

A{x,t + t), 

(7) 

A{-x,t), 

(8) 


for any (0, a, r) G G. In the present study it is the group G generated by the continuous 
symmetries (5)-(7) which enters the problem formulation. Namely, for a given solution 
A(x, t) of the CGLE, let us consider the isotropy subgroup Ga of G at A, 

Ga = {(<^, S,T)gG 1 Aix, t) = e^^A{x + S,t + T)}, (9) 

which consists of elements of the symmetry group G = x M leaving A invariant. With 
that in mind, we pose the problem: seek solutions A(x, t) of the CGLE satisfying 

A{x,t) = + S,t + T), (10) 

for (if, S,T) € G also unknown and to be determined. In other words, find orbits G • A 
of G generated by solutions A of the CGLE which are invariant under the action of some 
subgroup Ga C G, that is, Ga ■ A = A. Here, at least one subgroup of Ga generated by an 
element {ip, S,T) £ G is also to be determined. 
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As is clear from (10), the case (p = S = 0, T > 0, would result in a time-periodic solution. 
Within the more general context of the problem of seeking solutions of a dynamical system 
fixed by an action of its symmetry group which also contains time translation, as is the case 
resulting from T > 0 and nonzero or S' in (10), such invariant solutions are also referred 
to as relative time-periodic solutions. Since the solutions sought must satisfy the boundary 
(space-periodicity) condition (2), it is easy to see that if S = L^/q for some integer q> 1, 
then \A{x,t)\ = \A{x,t + qT)\, whereas if both S = L^/q and ip = 27r/g for some integer 
g > 1, then A{x,t) = A{x,t + qT) and, therefore, {0,0, qT) G Ga (i.e., A is time-periodic, 
with time period qT). 

Notice that if {ip,S,T) G Ga, the triples {jip,jS,jT), j G Z, are also elements of the 
isotropy subgroup Ga- Hence, {ip,S,T) generates a subgroup of Ga- Thus, the problem 
that we aim to solve numerically can be described succinctly as follows: 

1. Given a point po = {Ro,i'o, IJ-o) in the parameter space of the OGLE, find a solu¬ 

tion Apg{x,t) of the GGLE and a generator {ip{po), S{po),T{po)) of a subgroup of the 
isotropy subgroup such that condition (10) holds. That is, Ap^ is an invari¬ 

ant solution of the GGLE under the action of the subgroup of Ga^^^ generated by 
{ip{po),S{po),T{po)). 

2. Then, starting from p^ = {Rq,Vq, po), vary the point p = {R,v,p) along a subspace 
in the parameter space of the GGLE, ending at a point p^ = {R^,v^, p^), to find 
a sequence of new invariant solutions Ap{x,t) and generators {(p{p), S{p),T{p)) of 
subgroups of their corresponding isotropy subgroups Ga^,- 

In reference [35] we found 77 distinct unstable invariant solutions (that is, 77 G-orbits 
generated by distinct invariant solutions) of the GGLE at the point po = {Ro,Vq,Pq) = 
(16, —7, 5) of the parameter space of the GGLE, thus addressing the first part of the problem. 
Here, we take the first 15 of these solutions, per the listing from the Appendix in [35], and 
address the second part of the problem. Specihcally, using numerical continuation (as 
described in Section 3) we found 15 sequences (or discrete continuation paths) 

; (</?(pi*'),5’(pi*^),T(pi'^)) 1 

(11) 

pf = G [9,100] X [-7,-2.7] x [-0.05,5.98], 0 < A: < N«}, 

i = 1,..., 15, of new invariant solutions (10) of the GGLE and corresponding generators of 
subgroups of their isotropy subgroups Ga , . In (11), the number of invariant solutions 

in a sequence is at least 100 and, for each i = 1,..., 15, the final point in the GGLE 
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parameter space was fixed at ) = (100,-7,5). Thus, the sequences 

(11) can be thought of as a deformation of an initial set of distinct G-orbits at po = 
(16,—7, 5) into a final set of G-orbits at Pn = (100,-7,5), which in this study are also 
distinct with the only exception being that the final orbits in the sequences and A^*^ 
at Pn happened to coincide (details are provided in Section 4). In other words, if we 
think of the space of G-orbits as fibered over the parameter space of the CGLE (the base 
of the fibration), then the sequences (or continuation paths) A^’’'^ can be thought of as 
(discrete) sections of the fibration. Interestingly, several of the sequences A'-^^ that we have 
computed contain solutions with additional symmetries (which we describe in detail later 
in this section), thus revealing an intricate structure of the space of G-orbits, as well as 
corresponding bifurcation points on the base. 

Note that the meaning of the space-periodicity boundary condition (2) is that any 
solution A(x, t) in the class of solutions of the CGLE that we seek has a subgroup in its 
isotropy subgroup Ga which is generated by (0, L,,,, 0). In other words we restrict, a priori, 
the class of solutions of the CGLE that we look for to the one that contains, at a minimum, 
solutions with symmetry (2). This allows us to represent A(x,t) as a Eourier series 

A{x,t) = ^ a^(f)e‘*'™^, (12) 

mGZ 

where = 27rm/Lj, denotes the m-th wavenumber in the expansion. Erom the group- 
invariance condition (10) it then follows that the complex-valued Eourier coefficient func¬ 
tions a^{t) in (12) satisfy 

a^{t) = -h T) (13) 

for all m G Z. Because of the presence of symmetry (2), the solutions sought can be 
restricted to those with elements ((p, S,T) £ G having S G [0, L^). 

Moreover, since the CGLE is invariant under the action of the group Z 2 of spatial 
reflections A{x,t) A{—x,t), to any solution A{x,t) of the CGLE having ( 0 ,L 2 ., 0 ) and 
{ip,S,T) as generators of subgroups of the isotropy subgroup Ga (defined in (9)) there 
corresponds a solution A{x,t) ■= A{—x,t) having (0,L^,0) and {ip, L^—S,T) as generators 
of subgroups of the isotropy subgroup Ga- This can be seen from the chain of equalities 

A{x,t) := A{—x,t) = A{—x + S,t+ T) by (10) 

= e^'^A{x — S,t + T) by definition of A 

= e^A{x + {L,,-S),t + T) by (2). (14) 

To express the above in a more symmetric form, let us introduce 5 = \L,^/2 — Sj. Then, 
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if A{x,t) is a solution of the CGLE having (0,L^,0) and ± 6,T) as generators of 

subgroups of the isotropy subgroup Ga, the solution A{x,t) ■= A{—x,t) has ( 0 ,L 3 ,, 0 ) and 
{ip, L,^l2 =F 6, T) as generators of subgroups of the isotropy subgroup Ga- We shall call the 
invariant solutions {A-, ip, -L^/2 ± 5, T) and {A-, ip, L^/2 =F 5, T), as well as their corresponding 
orbits G ■ A and G • A, conjugate to each other under the (involutive) action of the group 
Z 2 of spatial reflection symmetry of the CGLE. 

Now, while invariance of solutions of the CGLE other than that defined by (10) and 
(2) is not part of the problem formulation (9)-(10), it is clearly not excluded from it. The 
CGLE may admit solutions having symmetries other than (or in addition to) that defined 
by (10) and several of the solutions resulting from our study do have additional symmetries. 
Thus, in what follows we discuss some such symmetries and their properties. We emphasize 
that our treatment on additional types of symmetries exhibited by solutions of the GGLE 
is not exhaustive, but rather inclusive of material relevant to the discussion on our results 
in Section 4. 

For instance, there may exist solutions of the CGLE satisfying 


A{x, 

.t) = 

e'^^/^A(x + L,,/l,t), 

for some 1 £ N, 1 > 1, 

(15) 

A{x, 

.t) = 

A{—x + 2ci, f) 

for some Ci G M, 

(16) 

A{x, 

.t) = 

-A{-x + 2c2,t) 

for some C 2 G M. 

(17) 


Note that (15) describes solutions fixed by a composition of the actions (6) and (5), and 
gives (2'k jl, L,^/l, 0) as one generator of a subgroup of Ga- From (15) it is also clear that the 
absolute value of such a solution has spatial period of L,^/l. Furthermore, by substituting 
condition (15) into the Fourier series expansion (12) it follows that the Fourier coefficient 
functions a™(t) of a solution with symmetry (15) satisfy 


{i) 


nonzero if m G {Irh — 1, m G Z} 
0 otherwise. 


(18) 


Symmetries (16) and (17) describe solutions that are, respectively, even about x = Ci or 
odd about x = C 2 for some real numbers Ci, C 2 . (These solutions are hxed by a composition 
of the actions (8), (6), and (5).) From (16) and the periodic boundary condition (2) it 
follows that a solution even about x = Ci is also even about x = Ci + Lj;/2; similarly 
a solution odd about x = C 2 is also odd about x = C 2 + 2. The Fourier coefficient 

functions a^{t) in (12) of a solution even about x = Ci satisfy 

a_^(f) = a^(f)e'''™^''i , m = 0,1,2,..., (19) 
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whereas for a solution odd about x = C 2 one has 

a_^(t) =m = 0,1,2,.... (20) 

From (18), (19), and (20) it follows that restricting the search for solutions to those possess¬ 
ing symmetries (15), (16), or (17) would lead to a reduction in the number of unknowns. 
However, as already mentioned, we did not make a priori such a restriction in order to 
allow for a more general set of solutions. Finally, we point out that a solution having both 
symmetries (15) and (16) also satisfies 

A(-x + 2(ci + Lj{2l)),t) = e‘2^/'H(x, t). (21) 

In particular, note that a solution satisfying (15) for I = 2 and which is even about x = Ci 
is also odd about x = C 2 = Ci -|- L^/A. 

For a solution satisfying (10) and (15) it follows that {ip — 2 'k/I, S — L^/l, T) is another 
generator of a subgroup of Ga- In particular, if it happens that for such a solution one 
has S = L^/l, then \A{x,t)\ = \A{x,t -|- r)| is satisfied, whereas if both S = L^jl and 
ip = 2tt/1, then A{x,t) = A{x,t + T) also holds. As for invariant solutions (10) that also 
possess symmetry (16), note that, for each m = 0,1, 2,..., 

a_„(t) = + T) by (13) 

= by (19). (22) 


On the other hand, for each m = 0,1,2,..., 

a.m.{t) = a™(t)e'''™^''i by (19) 

= e'V^"*-^a,„(t + r)e‘^’"2ci by (13). (23) 

From (22) and (23) it follows that we must have for all m G Z, which 

holds whenever S is an integer multiple of L„./2 (recall that = 21 ^ 1 /L^). The case 
for invariant solutions (10) with the additional symmetry (17) is analogous. Therefore, 
solutions satisfying (10) which also posses symmetries (16) or (17) exist in subspaces of 
the solution space {A] ip, S,T) for which either 5 = 0 or 5 = L^/2 (since, by the periodic 
boundary conditions (2), S can be restricted to be in the interval [0,L^)). 

To illustrate some of the aforementioned symmetries of solutions of the CGLE, Fig¬ 
ures 1-2 display several plots that aid in visualizing the invariant properties. Figure 1 
shows a solution of the CGLE having symmetry (10), but none of (15)-(17). This solution 
belongs to the sequence (see (11)) resulting from the numerical continuation procedure 
to be described in Section 3, that is, from the continuation path for the sequence listed 
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SR(A) 


«(A) 


K(A) 


!R(A) 


(a) f = 0 


(b) t = r/3 (c) t = 2r/3 


id) t = T 


X 

(e) 1^1 

Figure 1: Solution having symmetry (10), with T > 0 and nonzero and S. 

with id 5 in Tables 1-2. The time evolution, represented as curves on the plane with co¬ 
ordinates defined by the real part and imaginary part of the solution A{x,t) 

at different times within the interval [0, T], is depicted in Figures la-ld. The excitation 
of multiple temporal frequencies is apparent from these curves. For single-frequency so¬ 
lutions A{x,t) = B{x)e^‘^^ or generalized traveling waves A{x,t) = p{x — 

(where a; is some single frequency), plots of this kind would show, except for a rotation, 
the same curve at each point in time. Therefore it is clear that the solution depicted in 
Figures la-ld is not of either of these single-frequency types. Note also that the curve at 
time t = T differs only by a rotation from that at time t = 0 due to the rotation of the 
complex field A{x, 0) —>■ e^‘^A{x + S,T). Repeated patterns resulting from invariance due to 
time periodicity and the nonzero space translation S is better observed from surface plots 
of the absolute value |^| of A{x,t) over several space and time periods, as in Figure le. 

A solution possessing all of the symmetries (15)-(17), in addition to the symmetry 
(10), is shown in Figure 2. The plots represent a solution which belongs to the sequence 
of solutions under id 15 in Tables 1-3, computed at the point (R, i/,/x) = (100,-7,5) of 
the CGLE parameter space. Surface plots of the real part 3?(A), imaginary part A(A), 
and absolute value |A| of the solution A{x,t) are shown in Figures 2a-2c, where {x,t) G 
[0,2L^] X [0, 2T], that is, the surfaces are plotted over two space and two time periods. 
For this solution, symmetry (15) holds with I = 2. Since, in addition, the solution is even 
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Figure 2: Solution having symmetries (15), for I = 2, (16), and (17), in addition to (10). 


about X = mLj,/4, for m G Z odd, it follows from (21) that the solution is also odd about 
X = {rh + l)L,,,/4. Finally, the absolute value of the solution has spatial period L^l2 and is 
time-periodic, with period T. As seen in Figures 2a-2d, pattern similarities in both space 
and time are easily observed in the presence of the additional symmetries (15)-(17). 

We conclnde this section by noting the following fact. Suppose that for some {ip, c, T) G 
G, where G is the group of continuous symmetries of the CGLE (refer to (4)), a solution 
A(x, t) of the CGLE has the symmetry 

A{x,t) =e^^A{-x + c,t + f). (24) 

The Fourier coefficient functions a^(t) in (12) of such a solution satisfy 

a_m{t) = a^{t+ , m = 0,1,2,.... (25) 

The right-hand side of (24) is the result of the (left) action on A(x, t) of the composition 
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{(p,c,T) o (x —)■ —x), and after a (left) action of said composition on both sides of (24) one 
obtains that 

A{x,t) = e^'^A{—x + c,t + T) 

= e‘2‘^A(x,t + 2r). (26) 

Hence (2(^,0, 2T) G Ga, where Ga is the isotropy subgroup defined in (9). Also, note that 
we have 

[{(p,c,f) o {x ^-x)f = (2(^,0, 2f) G G, 

for any element {(f, c, T) in the group G of continuous symmetries of the CGLE. Conversely, 
let {(fi, 0, T) G Ga for some solution A{x, t) of the CGLE. Then we have 

[((^/2 + fcvr, c,T/2) o (x — )•-x)]^ = 

for every k € Z and c G M. Therefore, solutions of the CGLE with symmetry (24) do 
possess symmetry (10) of the type we seek, and, conversely, a solution with symmetry 
(10) may also possess the additional symmetry (24). An example of solutions having both 
symmetries (10) and (24) is described in Section 4 (cf. Figure 9). Such solutions appeared 
in the continuation path for the sequence listed with id 13 in Tables 1-3, that is, in the 
sequence (see (11)). 

Finally, note that if a solution A{x,t) of the CGLE having symmetry (10) for S' = 0 or 
S = L^/2 also satisfies 

A(x, t) = e^‘^A(—X + c, t) (27) 

for some real numbers ip and c (an example being solutions with the additional symmetry 
(16), where (^ = 0, or with the additional symmetry (17), for which <p = tt), then 

A{x,t) = e‘(‘^+2c^)A(-x + c,t + 2 r). 

That is, such solution A also has symmetry (24). Therefore, one should expect to find 
solutions having both symmetries (10) and (24) in subspaces of the space of solutions 
{A-,ip,S,T) for which either 5 = 0, by (26), or 5 G {0, L,,./2}, if symmetry (27) is also 
present. 

3 Numerical Method 

As noted in Section 1, having computed previously in [35] a set of unstable invariant 
solutions of the GGLE for fixed values of the parameters {R, u, //), our first goal is to employ 
numerical continuation to carry solutions of this initial set into solutions in a regime with 
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a different set of parameter values {R,i',n). The methodology used towards this end is 
explained next. 

3.1 Discretization 

We employ Fourier series expansions in both space and time to derive an underdetermined 
system of nonlinear algebraic equations from which invariant solutions of the CGLE are 
sought. Sections 3.1.1-3.1.2 discuss the derivation of this underdetermined system, some 
of its properties, and those of its Jacobian matrix. 


3.1.1 Derivation of Nonlinear Algebraic Equations 


Since the boundary conditions (2) are periodic in x, we use the spatial Fourier series (12) and 
substitute into the CGLE (1) to obtain an infinite system of ordinary differential equations 
(ODEs), 

-^ = Ra^ - - {I + ifj.) ^ (28) 

miH-m2—m3=m 

for the complex-valued functions a^{t). Under this transformation the symmetries (5)-(8) 
of equations (l)-(2) become symmetries of (28). Thus, if a{t) = {a^{t)) is a solution of the 
system of ODEs (28), then so are 


(e‘^a„(t)). 

(29) 

V™'^a„(t)), 

(30) 


(31) 

(^ —m(^))? 

(32) 


for any {9, a, r) G x M. In particular, (29) and (30) say that the ODEs (28) are invariant 
under the T^-action 

{9, a) ■ (a„(t)) = (e‘®e'™‘^o,„(t)). 

In our computations we employ a spectral-Galerkin projection obtained by fixing an 
even number Nx and truncating the expansion (12) to include only the terms with indices 
m satisfying —Nxl2 + 1 < m < Nx/2 — 1. We then work with the corresponding finite 
system of ODEs which results from (28) after the Galerkin projection. Much accumulated 
theory and computation shows that for sufficiently large Nx the behavior of this truncation 
captures the essential features of the dynamics of (l)-(2) [16, 25]. The system of Nx — l 
complex ODEs resulting from the Galerkin projection has the same symmetries (29)-(32) 
as the infinite system. The T^-action is linear on the finite-dimensional space with 
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a group element {9, a) G acting via multiplication by the unitary matrix diag(e‘®e^™''^), 
where it is understood that the indexing of the matrix elements is done by increasing value 
of m = —Nx /2 +1,..., Nx /2 — 1. Henceforth this indexing is assumed whenever applicable. 

From the condition (10) defining an invariant solution of the CGLE, it follows that the 
corresponding solution a(t) of the system of ODEs (28) satisfies 

a^(f) = e“^e‘^™'^a„(f + T) (33) 

for all m and t (and where ip^ S,T are to be determined). It is easy to see that the set of 
functions 

aUt) = ^ , (34) 

riGZ 

where = 2TTnjT denotes the n-th frequency in the expansion, are a solution of the 
system of functional equations (33). Hence, they provide an appropriate representation for 
invariant solutions of the system of ODEs (28). Substituting (34) into the truncated system 
of ODEs (28) and using again a Galerkin projection obtained by fixing an even number Nt, 
so that the summation index in (34) runs over the range —Nt/2 + 1 < n < Nt/2 — 1, results 
in a system of nonlinear algebraic equations. 



T 



- A;^(l + ii^)a,„,„ 


(1 + i/r) 


E ( E 

mi+m2—rrL^=rrL \ni-\-n2—nz=n 


a 


m-^ ,n-^ 


a 


7712,n2 



(35) 


for the complex Fourier coefficients and elements {(p, S, T) of the isotropy subgroup 

(9). Note that (35) is an underdetermined system of {Nx — ^){Nt — 1) complex equations in 
(Na; —l)(Nt —1) complex unknowns plus three real unknowns or, after splitting the equations 
into their real and imaginary parts, 2[Nx — ^){Nt — 1) real equations in 2{Nx — ^){Nt — 1) + 3 
real unknowns. Solutions of this system of equations will give the desired invariant solutions 
of the truncated system of ODEs via the expansion (34). We note here that with the 
introduction of the representation (34), the symmetry group G = x M of (1) and (28) 
descends to the symmetry group x of (35), acting on the space ({a^^„}, p, S, T) 

of solutions of (35). Henceforth, by a slight abuse of notation, we refer to both symmetry 
groups X M and as G. 

For convenience, we separate (35) into its linear and nonlinear parts and write it as 

F(a, (^,5',T) = FL(a,(^,5',r)+ FNL(a) = 0, (36) 
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where d is a vector with components given by the coefficients the vector FL(d, ip, S, T) 

has components given by 

(^i ^ -R + ki{l + a^,„, 

and FNL(d) is a vector with components given by 

(1 + i^) E E ^7712 ,Ti-2 ^TTT-3 l * 

mi+m2—rM3=m \ni+n2—n3=n / 

Note that the components of FNL(d) are simply the coefficients in the truncated Fourier 
series expansion (in both space and time) of the function (1 + ifi)A\A\‘^. Furthermore, in 
defining the vector a (and similarly for Fl and Fnl) we are implicitly assigning an ordering 
on the coefficients {d,„,„} that uniquely determines an indexing for the components of a. 
Henceforth, such a convention should be understood whenever applicable. Finally, we will 
use the notation in (36) to denote both the system of complex equations (35) and the 
system obtained by splitting (35) into its real and imaginary parts, as it should be clear 
from the context which case applies. 

The symmetries (29)-(32) of the system of ODEs (28) induce symmetries of the system 
of nonlinear algebraic equations (35). Note that if S,T) is a solution of F = 0, 

then for any {6, a, r) G 


({e‘'^ar„,„},(/?, 5, T), 

(37) 

({e'™"d„,„},(^,5,r). 

(38) 


(39) 

S', T), 

(40) 


are also solutions. From the continuous symmetries (37)-(39), it follows that the set of 
solutions of F = 0 splits into orbits st) symmetry group T^, 

^(a,ip,s,T) ■= • {a,ip,S,T) \ (0,a,T) € T^} , (41) 

where the action of on a point (d, p, S, T) is defined by 

{e,a,T)-[a,p,S,T) = ({e'V™V""a„,,J,v^,5,r). (42) 

That is, acts on d via multiplication by the matrix diag(e'^e™'^e'”'’'), and it acts trivially 
on (y?. S', T). Finally we note that, for the system of nonlinear algebraic equations (35), the 
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transformation (40) induced by (8) maps a solution 

{{a^,n},‘f,L^/2±5,T) (43) 

of (35) to another (conjugate) solution 

of (35), where, again, 5 = \L^/2 — S'], (Refer to the paragraph containing (14).) We now 
proceed with a brief discussion concerning the Jacobian matrix of the system (35). 

3.1.2 Jacobian Matrix 

One can see from the system of equations (36) that, although the individual components 
of the Jacobian matrix of F coming from the discretized linear terms FL(d, 5, T) in the 
CGLE are easily computed, such is not the case for those arising from the discretization 
of the nonlinear terms FNL(d). Hence explicit computation of the full Jacobian matrix 
from (35) would be cumbersome. In addition, the Jacobian matrix is dense, so, as the 
number of unknowns (and equations) increases, it becomes unfeasible to solve linear systems 
with the Jacobian as coefficient matrix using direct methods. However, matrix-vector 
products with the Jacobian matrix of F can be computed efficiently for the problem at 
hand, making the use of iterative methods for solving linear systems a reasonable option. 
We proceed to review the calculation of this matrix-vector product since, as will be discussed 
in Section 3.2.1, it is an essential feature of the Newton step computation employed in the 
numerical continuation. 

Let denote the matrix whose columns correspond to derivatives of F with respect to 
the real and imaginary parts of the unknowns and let h be a vector with components 

given by the coefficients in the truncated Fourier series expansion of a function V{x,t). 
Assume that J- is evaluated at a given point {a, (p, S,T). The product J^v can then be 
computed as^ 

J-u = DFL(h,y5, S'jT)-4 DFNL(a,'y), (45) 

where DFl(u, ip, S, T) = FL(f), (/?, 5, T) (as defined in (36)) and DFNL(d, v) is a vector with 
components given by the coefficients in the truncated Fourier series expansion (in both space 
and time) of (1 -|- ipL){A^V* + 2\A\^V). This matrix-vector product operation follows from 
the discretization (analogous to that used for the CGLE) of the first variational derivative 

^Note that in the right-hand side of (45) we are actually using v to denote a vector with the complex 
numbers {vm,n} as components, whereas in the left-hand side of (45) v denotes a vector with real components 
that are the real and imaginary parts of the coefficients {Vm,n}- We make use of this slight abuse of notation 
in this paper since the intended meaning should be clear from the context. 
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of equation (1) 


(l + z/i)(AV* + 2|^|V). 


dV d^V 

^ = w + (i + .o^- 

Furthermore, as can be seen from system (35), the operation of computing a matrix-vector 
product with the columns of the Jacobian matrix of F corresponding to the derivatives 
with respect to (/?, S', and T poses no difficulty. 

It follows then that matrix-vector products with the Jacobian matrix of F can be easily 
computed without the need of explicitly calculating the (full) Jacobian. Note also from 
(36) that the portion of the Jacobian matrix coming from the discretized linear terms 
FL(d, (/3, S, T) in the CGLE is a block diagonal matrix, with 2x2 blocks, whose components 
are easily computed. Hence, solving linear systems with this block diagonal matrix poses no 
complications. This is advantageous since, as will be discussed in Section 3.2.1, this block 
diagonal matrix provides an effective preconditioner for some of the matrix-free iterative 
methods [23, 28] when solving linear systems having the Jacobian as coefficient matrix for 
the problem at hand. 

Now recall that, split into its real and imaginary parts, the system F = 0 is an underde¬ 
termined system with three more unknowns than the number of equations. Thus, the kernel 
of the Jacobian matrix of F, ker(J), is at least three-dimensional. From the T^-invariance 
of F (due to the symmetries (37)-(39)), one obtains three linearly independent vectors in 
ker(J) at a solution {a,if,S,T) of F = 0. Such vectors result from a basis for the space 
of infinitesimal generators of the action (42) of on the point (d, (^, S, T), which can be 
obtained as 

d 

d« 

( 


*3 := ^ ({e- 



= ({id„,„},0,0,0), 

J 

0=0 


= ({imd,„,„}, 0,0,0) 


(T=0 


= ({in d„,,„}, 0,0,0) 


r=0 


Split into real and imaginary parts, the above vectors vi, V 2 , and are in ker(J) at a 
solution (d, </?, S, T) of F = 0. Note, in particular, that the matrix J- (refer to (45)), whose 
columns correspond to derivatives of F with respect to the real and imaginary parts of the 
unknowns {d™,„}, is singular at a solution of F = 0. We will come back to this point in 
Section 3.2.1, where the computation of the Newton step is addressed. This concludes the 
discussion on the discretization; details related to the numerical continuation now follow. 
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3.2 Numerical Continuation of Solutions 


The numerical continuation was done using the Library of Continuation Algorithms (LOCA) 
software package [46], specifically with the aid of the algorithms provided to track steady 
state solutions of discretized PDEs as a function of a single parameter. The option of 
arc-length continuation was used in order to allow for turning points in the path following 
process. Although we are not computing steady state solutions in this study, it is clear 
that the feature of tracking steady state solutions in the LOCA package provides the ca¬ 
pability of solving a system of nonlinear algebraic equations using numerical continuation. 
We thus take advantage of this feature, particularly to handle the step size control in the 
continuation parameter. 

Newton’s method is used to solve the system of nonlinear algebraic equations in the 
numerical continuation, and the user must supply the LOCA package with a routine for 
computing the Newton step. That is, the user must provide a routine that solves a linear 
system having as coefficient matrix the Jacobian of the system of nonlinear algebraic equa¬ 
tions. For this purpose, we employed iterative methods for solving linear systems [23, 28], 
specifically solvers for linear systems from the Meschach software package [47]. For fur¬ 
ther efficiency in the calculations, we used POSIX threads (pthreads) programming [36] in 
our routines, taking thus advantage of the multiple cores available nowadays in personal 
workstations. Relevant details about this Newton step computation are provided next, in 
Section 3.2.1. Finally, the computation of the nonlinear terms Fnl in (36) and DF^l in 
(45), needed, respectively, for the evaluation of F and that of the product of the Jacobian 
matrix and a vector, was done using the FFTW software package [20]. 

3.2.1 Newton Step Computation 

Rather than augmenting the system (35) with an additional set of equations in order to 
work with an equal number of equations and unknowns [35], we work with the underde¬ 
termined system (35) and consider here a Newton step defined from the Moore-Penrose 
inverse [6], which yields a minimum norm solution of the system of linear equations with 
the underdetermined Jacobian as coefficient matrix, and is one technique used in numeri¬ 
cal continuation methods [1, 54]. However, instead of computing the desired Newton step 
directly from the linear system having the underdetermined Jacobian as coefficient matrix, 
we premultiply the linear system with a matrix composed of a subset of the columns of the 
Jacobian so that a numerical solution for the problem at hand may be obtained in a more 
efficient manner, as explained next. 

Let J be a p X g matrix, p < q, and assume J has rank p. Let 6 be a vector of size g x 1. 
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Recall that the minimum norm solution z of the system of linear equations^ 

Jz = b (46) 

given by the Moore-Penrose inverse is [1] 

z = (47) 

Computing the solution z from (47) thus requires the solution of a linear system having 
as coefficient matrix. As is well known, the numerical solution of a system of linear algebraic 
equations may be obtained using a variety of methods, either direct [21] or iterative ones 
[23, 28]. A main distinction between these two classes of methods is that the use of direct 
methods requires explicitly the availability of the coefficient matrix, whereas for iterative 
methods what is needed is the ability to perform matrix-vector products with the coefficient 
matrix. Hence, iterative methods are an attractive option when explicit computation of 
the coefficient matrix is not feasible or is inconvenient, and multiplication of the coefficient 
matrix and a vector can be performed (efficiently) without explicit computation of the 
coefficient matrix. As discussed in Section 3.1.2, the latter applies to the problem at 
hand. Therefore we considered the use of iterative methods, in particular the generalized 
minimal residual (GMRES) method [23, 28, 45] due to its robustness and suitability for 
non-symmetric systems. Since the matrix in (47) is symmetric we also explored the 
possibility of using the conjugate gradient (CG) method for symmetric systems [23, 28], 
but, as will be discussed below, the GMRES method is a more suitable option for our 
problem. 

The convergence behavior of iterative methods for solving linear systems is dependent 
on the method as well as on various other factors, for example, certain properties of the 
coefficient matrix or the problem from which the linear system is derived. In the case of 
the GMRES method, one desirable property for fast convergence is for the eigenvalues of 
the coefficient matrix to be clustered around a few values, away from zero [23, 28]. Another 
important issue is that the use of iterative methods for solving linear systems typically 
requires the use of a preconditioner in order to perform efficiently. Generally speaking, 
preconditioning refers to multiplying the linear system (on the left or right) by a matrix such 
that the resulting system has the properties needed for optimal or enhanced performance of 
the particular method under consideration (and yields a solution for the unpreconditioned 
(i.e., original) linear system). For thorough treatments on iterative methods for solving 

^The components of the Jacobian matrix in our computations are real, hence our use of the terms 
transpose and symmetric when referring to the matrix J in the discussion that follows, instead of the more 
general terms conjugate transpose and Hermitian. 
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Figure 3: Typical convergence behavior of GMRES and CG for the problem at hand. 
Suitable convergence behavior results for GMRES with DgJs as coefficient matrix, and 
for GMRES or GG with as coefficient matrix. However, the only viable option 

for our problem is the GMRES method with as coefficient matrix. 


linear systems, the interested reader is referred to [23, 28] and references therein. Here we 
restrict ourselves to a brief discussion on the behavior resulting from the use of the GMRES 
and GG methods for the problem at hand. 

Eigure 3 depicts the typical convergence behavior exhibited by the GMRES and CG 
methods when used to compute the solution of systems of linear equations having as coeffi¬ 
cient matrix the Jacobian of the system (35) of nonlinear algebraic equations, after splitting 
the equations into their real and imaginary parts. In the figure, J denotes the underde¬ 
termined Jacobian matrix of the system (35), Jg represents a square, non-singular matrix 
composed of a subset of columns of J, the preconditioner Dg is the Jacobian matrix of 
the discretized linear terms Fl (36), with columns corresponding to the derivatives with 
respect to the real and imaginary parts of the unknowns (it is a block diagonal 

matrix), and the preconditioner U is a diagonal matrix having the diagonal elements of 
JJ^ on its diagonal. As can be seen from Figure 3, suitable convergence behavior results 
only from the use of GMRES with Dg^Jg as coefficient matrix (i.e., the use of GMRES to 
solve linear systems with Jg as coefficient matrix and Dg as a preconditioner), as well as 
from using GMRES or CG with D~^ JD~^ as coefficient matrix (i.e., the use of GMRES 
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or CG to solve linear systems with JJ"^ as coefficient matrix and D as a preconditioner). 

The block diagonal preconditioner Dg is very effective when used with the GMRES 
method to solve linear systems with Jg as coefficient matrix, as seen in Figure 3. For 
the example depicted, it took 190 iterations to solve for a system having 5,925 unknowns. 
The corresponding run with unpreconditioned GMRFS required 5,676 iterations, making 
it impractical for our purposes. Using the preconditioner D and both the GMRFS and CG 
methods to solve for systems having JJ^ as coefficient matrix also gave good results. For 
the example in Figure 3, these two methods took, respectively, 269 and 552 iterations to 
solve the preconditioned system. However, per the discussion in Section 3.1.2, the block 
diagonal preconditioner Dg is readily available and easy to manipulate, whereas assembling 
the preconditioner D requires calculating the diagonal terms of the matrix JJ^, and these 
terms are not readily available for our problem. Therefore, the only viable option for us is 
the use of the GMRES method with preconditioner Dg to solve systems having Jg as the 
coefficient matrix. As a result, the computation of the minimum norm solution z directly 
from (46)-(47) is unfeasible for the problem at hand. 

Returning to the linear system in (46), we thus express J as composed by two matrices, 

J=[Jg\Jr], (48) 

where Jg has dimension p x p (and is non-singular) and Jr has dimension p x {q—p). Now 
we consider the system z = obtained by multiplying (46) with yielding 

[l\J-^Jr]z = J-\ (49) 

where I is the p x p identity matrix. We work directly with the system (49) and compute 
the desired solution z by solving two sub-problems, namely: 

1. Compute the right-hand side Jj^b, as well as the q—p columns of the submatrix Jj^Jr 

of the coefficient matrix [/| in (49). This sub-problem will therefore require 

the solution of q—p-\-l linear systems having Jg as coefficient matrix. It will be feasible 
if solving linear systems having Jg as coefficient matrix can be done efficiently and if 
q—p+1 is small. Both of these requirements are satisfied in our study since first, per 
the discussion from the previous paragraphs, we can use the GMRES method to solve 
the required linear systems efficiently, and second, for our problem, q—p = 3 due to 
the 3-tuple {(p,S,T) of additional unknowns in the problem formulation (9)-(10). 

2. Upon completion of sub-problem 1, compute the minimum norm solution given by 
the Moore-Penrose inverse for the underdetermined system of linear equations (49) . 


23 


Per sub-problem 1 above, one first needs to solve A: = 1,..., q—p+1 linear systems 

JsZk = bk, (50) 

where, for k = 1,..., q—p, the A:-th linear system (50) will have the right-hand side vector 
equal to the A:-th column of the matrix Jr, so that the solutions i*. of said q — p linear 
systems yield the columns of the submatrix Jj^Jr of the coefficient matrix [/| in 

(49). The solution of the remaining linear system (50), with bq_^i = b, yields the right-hand 
side Jj^b in (49). Solving these q—p+1 linear systems (50) is not an obstacle since they 
can be solved either independently, in parallel, or with an implementation of a (direct or 
iterative) solver for linear systems that handles multiple right-hand sides. Recall that the 
viable option for us is to use the GMRES method for linear systems. In our implementation, 
we combined it with the use of POSIX threads (pthreads) programming [36] in order to 
solve the required linear systems (50) in parallel. Hence, in the current study, the solution 
of the q—p+1 = 4 linear systems (50) having Jg as coefficient matrix was achieved basically 
in the same amount of time as that required to solve a single such system. 

Upon completion of sub-problem 1, what remains to be done is to compute the minimum 
norm solution z from the system in (49). As previously noted, one desirable property for 
fast convergence of the GMRES method is for the eigenvalues of the coefficient matrix to 
be clustered around a few values, away from zero, since, typically, the number of iterations 
required for convergence when using the GMRES method depends on the number of distinct 
eigenvalues of the coefficient matrix of the linear system [23, 28]. Thus, we also used the 
GMRES method to solve for the minimum norm solution 2 : in (49), since this requires 
solving a linear system having 

I+{J-^Jr)iJ-^Jrf (51) 

as coefficient matrix, and the matrix (51) has all but q—p eigenvalues equal to one, with 
the remaining eigenvalues greater than or equal to one. (Recall that Jg is a non-singular 
px p matrix and Jr has dimension p x {q—p), where p < q.) This follows from the fact that 
the eigenvalues A and eigenvectors y of the matrix (51) satisfy 

{J-^Jr){J-^Jr)^y = (A-l)y. (52) 

Noting that the null space of the matrix {J+^Jr)^ has dimension (at least) p— {q—p) = ‘^p—q, 
it follows that the matrix {J+^ Jr){J+^ Jr)'^ has zero as an eigenvalue, that is, A = 1, with 
multiplicity (at least) 2p — q. Furthermore, {J+^Jr){J+^Jr)"^ is positive semi-definite and 
symmetric, so its eigenvalues (A —1) in (52) are non-negative and, thus, the remaining (at 
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most) q—p eigenvalues of the matrix (51) satisfy A > 1. The significance here is that solving 
a linear system with the matrix (51) as coefficient matrix, which is required in order to 
compute the minimum norm solution z of system (49), should take q—p+1 iterations if 
we use the GMRES method. For our problem, this means q—p+1 = 4 iterations. Note 
also that computing matrix-vector products with the matrix (51) can be easily done (since 
the q—p = 3 columns of the matrix Jj^Jr have been previously computed and stored in 
memory). Furthermore, no preconditioning is required to solve linear systems having (51) 
as coefficient matrix. Hence, using the GMRFS method to solve the aforementioned sub¬ 
problem 2 poses no difficulty and results in a negligible amount of additional computing 
time when solving for the Newton step using the proposed approach. 

Finally we note that the vector h in (49) corresponds to —F evaluated at the current 
solution estimate (from Newton’s method). Also, based on our presentation of the material, 
it may seem natural to consider the matrix J- introduced in Section 3.1.2, whose columns 
correspond to derivatives of F with respect to the real and imaginary parts of the unknowns 
{dm.n}, as that corresponding to the matrix Jg in (48)-(49). Recall from Section 3.1.2, 
though, that J- is singular at a solution of F = 0. We therefore define the matrix Jg 
as that obtained by replacing three columns^ from J- by the columns of the Jacobian 
matrix of F corresponding to derivatives with respect to the unknowns (p, S, and T. This 
proved effective in dealing with said singularity when computing the Newton step from 
the corresponding system (49) during the numerical continuation. The replaced columns 
from J' define the columns of the matrix Jr in (48), which are needed to construct the 
matrix Jj^Jr in (49). These three columns can be computed (efficiently) via matrix-vector 
products of the Jacobian matrix (see Section 3.1.2) and standard basis vectors. 

4 Numerical Study and Results 

The procedure described in Section 3 was applied to a subset of the unstable invariant 
solutions of the GGLF computed at the point {R, u, p) = (16, —7, 5) of the CGLF parameter 
space (without employing continuation) in the preceding study [35] , in order to carry them 
into solutions in a regime with an increased value of the parameter R, namely to the 
region at the point = (100,-7,5) of the GGLF parameter space. As indicated in 

Section 1, chaotic behavior is exhibited both at the initial and final parameter regions. 

A summary of the obtained results is gathered in Tables 1-3 and Figure 4b below. 
Already from them, we see that our probe into the moduli space of G-orbits reveals a 

"^As discussed in [3.5], the kernel of the Jacobian matrix at a solution of F = 0 is typically three- 
dimensional. 
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complicated and interesting structure. To start with, along each of the continuation paths 
A^''\ i = 1, , 15, we were able to compute a number (listed in Table 2) of new distinct 
G-orbits of solutions of the CGLE (each one of which corresponds to a distinct invariant 
solution of the CGLE). Thus, the continuation paths i = 1,... , 15, can be thought of 
as (discrete) sections of the fibered space of G-orbits over the space of parameters of the 
CGLE. 

Before describing the content of Tables 1-3, let us list the possibilities that may occur 

when numerically continuing a set of distinct G-orbits. (These are analogous to the cases 

listed later in Section 4.1 in which a point in the CGLE parameter space returns to itself 

after a series of steps while performing continuation for a single G-orbit.) Suppose that 

A (i) G and A q) G A^^\ i A Ji invariant solutions representing distinct G-orbits 

f’o Pq 

at the initial point po^ = Po ^ = (Ro, i^o, fJ-o) in the CGLE parameter space, where A^'’'’ and 
A^^^ are, respectively, the continuation paths (11) emanating from each one of the two initial 
invariant solutions. Given two points p^A and pA in the GGLE parameter space and two 
invariant solutions A (u G and A (p G A^^\ it may happen that p^A = pA , and we have 

Pk Pi 

to consider several possibilities. Namely, whether the invariant solutions A^(i) and 
represent (i) the same G-orbit, that is, {(fiPk^), S{pl^^),T{pi"'')) = {^{pY^), S{pA),T{pY'’)) 
and there exists some {9,a,T) G G such that A {e, = {0,a,T) ■ A(j)-, (ii) different, but 

Pk Pi 

conjugate, G-orbits, as defined in Section 2 (see also (43)-(44)); or (iii) different, non¬ 
conjugate, G-orbits. 

After a careful analysis of all computed solutions we found that there were solutions 
A^{ 2 ) G A^^'' and A^{4) G A^'^^ which represent the same G-orbit, at points pf^ = pY'' in 
the range {R,i^,p) G ([16.5,54.4] U [54.7,100]) x [—7] x [5] of GGLE parameter values. 
Furthermore, there were (a) solutions A (n) G and A (i4) G which represent 

Pk Pi 

the same G-orbit, at points pY^^ = pY*^ in the range {R,y,pL) G [11.5,60] x [—7] x [5]; 
(b) solutions A (ip G A (i4) G and A (is) G which represent the same 

Pk Pi Pj 

G-orbit, at points pY^^ = pY^^ = in the range {R,^,^) G [16,60] x [—7] x [5]; (c) 
solutions A (i4) G and A (is) G A^^®' which represent the same G-orbit, at points 

'Pj 

pY'^^ = pY^'’ in the range {R,u,p) G [16,85] x [—7] x [5]; and (d) solutions A (s) G A^®' 

Pk 

and A (6) G A^®' which represent conjugate G-orbits, at points p® = pf'’ in the range 

Pi 

{R,u,p) G [17,59.5] X [—7] x [5]. The latter case is illustrated in Figure 4a, where the 
graphs in (R, r)-space for the sequences A^®’ and A^®' overlap for the aforementioned range 
of R. Therefore, multiple representatives of same G-orbits were carefully accounted for and 
only one of them was taken as representative of the corresponding distinct G-orbit. 

Figure 4a illustrates as well that, while for the sequences A^^' and A^®\ for example, the 
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(a) continuation paths in {R, T) space 



continuation sequence id 

(b) unstable dimension 


Figure 4: (a) Representative continuation paths, depicted by plotting the time period T as 
a function of the parameter R. (b) Unstable dimension of solutions at the initial and final 
parameter regions for each sequence A^''\ i = 1,..., 15. 


numerical continuation progressed in a relatively smooth manner, such was not the case in 
general. Turning points and intersecting paths, exemplified by the depiction of the graphs 
T = T{R) for sequences and in Figure 4a, were frequently encountered. 

These features revealed intricate and challenging parameter regions for traversal. (Details 
appear in Section 4.1 below.) 

Table 1 lists the values of (ip, S, T) for the solutions at the starting CGLE parameter 
values of {R,v,r) = (16,—7, 5) and at the final values of {R,iy,fi) = (100,-7,5), as well 
the unstable dimension^ and spatial period of the invariant solutions at the aforementioned 
parameter values. Per the third column in Table 1, the listed solutions are unstable. As 
seen from Table 1 and the depiction in Figure 4b, the solutions used as initial points for 
the numerical continuation have unstable dimension ranging between 3 and 6, whereas the 
new solutions in the final parameter region have unstable dimension between 14 and 23. 
The time period for the initial solutions is in the range T G (0.02, 0.12) (or T G (0.32,1.92) 
for the formulation (3) of the CGLE); for the new solutions in the final parameter region 
we have T G (0.001,0.11) (or T G (0.1,11) for the formulation (3) of the CGLE). No truly 
time-periodic solutions were identified (although their existence in the regions traversed is 
not ruled out), as all solutions have a nonzero value for the rotation angle ip. 

®The unstable dimension of an invariant (or relative time-periodic) solution is the number of eigenvalues 
of the associated relative monodromy matrix having magnitude greater than one; see [35]. 
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id 

{ip, S, T)(^ 7 5) —S' {if, S, T)(^ ,, ,^) = (i00,-7,5) 

unstable 

dimension 


1 

(5.3622, 3.8544, 0.0233) —> (5.9158, 3.8856, 0.0015) 

4 —>16 


2 

(2.8849, 3.0956, 0.0539) —^ (0.1088, 3.1416, 0.0130) 

5 —>21 


3 

(0.0011, 3.9709, 0.0539) —^ (5.5905, 2.2876, 0.0193) 

5 —>22 

^x 

4 

(2.9343, 3.1416, 0.0540) —^ (0.1088, 3.1416, 0.0130) 

4 —>21 

Lx 

5 

(4.6093, 1.4537, 0.0547) —^ (0.9483, 1.0333, 0.0567) 

5 —>20 

Lx 

6 

(4.5165, 4.7061, 0.0556) —^ (5.2620, 5.1417, 0.0374) 

5 —>23 

Lx 

7 

(0.2436, 2.3887, 0.0608) —^ (4.0066, 3.1416, 0.0319) 

4 —>18 

Lx 

8 

(4.7959, 3.0824, 0.0825) —^ (3.8358, 3.4537, 0.0859) 

5 —>22 

Lx 

9 

(0.2876, 2.4431, 0.0875) —^ (0.3410, 1.3964, 0.0047) 

6 —>14 

Lx 

10 

(5.0251, 3.1416, 0.0895) —^ (5.3410, 3.1728, 0.0491) 

3 —>20 

Lx 

11 

(2.6023, 3.1719, 0.1046) —^ (1.5060, 3.2037, 0.0762) 

4 —>21 

Lx 

12 

(2.6575, 3.1209, 0.1078) —^ (4.5024, 2.4768, 0.0754) 

3 —>18 

Lx 

13 

(6.0553, 0.0032, 0.1106) —^ (2.5186, 0.0000, 0.0803) 

6 —>20 

Lx 

14 

(2.6063, 3.1057, 0.1128) —^ (4.0182, 3.2164, 0.0948) 

4 —>18 

Lx 

15 

(2.2500, 3.1416, 0.1146) —^ (1.7332, 3.1416, 0.1020) 

3 —>18 

Lx 


spatial 

period 





Table 1: Properties of solutions at initial and final points of continuation. 


Except for the sequence listed with id 1 in Tables 1-2, for which the solution at 
the final parameter values has only a few temporal frequencies active and appears to be 
close to a single-frequency solution, all of the solutions have broad spatial and temporal 
spectra. Also, aside from the sequences and listed, respectively, with ids 2 and 4 
in Tables 1-3, all of the resulting solutions retained the same spatial period of length 
as that of the starting solutions. The spatial period T^/3 of solutions in the sequences A^^'’ 
and A^'^^ was acquired (for both sequences) at parameter values ^ (20.2,—7, 5). 

The ending solutions in these two sequences are different elements of the same orbit (41) of 
the symmetry group G at the final point in parameter space, although the corresponding 
starting solutions belong to different orbits. As for the other sequences, the solutions at 
the hnal point in parameter space belong to different orbits of the symmetry group G. 

Breaking or gaining of the additional symmetries (16) or (17) was often detected, and 
gain of the additional symmetries (15) and (24) was also uncovered. (More details appear 
in Tables 2-3 and Section 4.1.) We did not observe a change in stability of the solutions 
at the points where the symmetry changes (the solutions remained unstable), but the 
unstable dimension would usually change at said points (with an increase or decrease of 1 
or 2). Table 2 indicates which additional symmetries, if any, the invariant solutions posses, 
whether continuation was done only on the parameter R or not (as will be discussed in 
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additional symmetries 

continuation 


id 

start of continuation 

in between 

end of continuation 

on R only 

N(') 

1 

none 

none 

none 

yes 

112 

2 

none 

(16) 

(16) 

yes 

188 

3 

(15), 1 = 2 

(15), 1 = 2 

(15), / = 2 

no 

179 

4 

(16) 

(16) 

(16) 

yes 

233 

5 

none 

none 

none 

yes 

634 

6 

none 

none 

none 

yes 

191 

7 

none 

(15) , / = 2, 

(16) , (17) 

(15) , / = 2, 

(16) , (17) 

no 

179 

8 

none 

(16) 

none 

yes 

489 

9 

(15), 1 = 3 

(15), 1 = 3 

(15), 1 = 3 

yes 

116 

10 

(16) 

(16) 

none 

no 

385 

11 

none 

(15) , / = 2, 

(16) , (17) 

(15), / = 2 

yes 

415 

12 

none 

none 

none 

no 

472 

13 

none 

(24) 

(24) 

no 

526 

14 

none 

(15) , / = 2, 

(16) , (17) 

(15), / = 2 

yes 

615 

15 

(15) , / = 2, 

(16) , (17) 

(15) , / = 2, 

(16) , (17) 

(15) , / = 2, 

(16) , (17) 

yes 

438 


Table 2: Types of additional symmetries associated to G-orbits along the continuation 
paths. 


Section 4.1), as well as the number of distinct CGLE parameter points for which 
solutions were found in each sequence i = 1,...,15 (see (11)). To determine the 
number we counted two points in the resulting numerical continuation path of the 
sequence , say p)*' = and = {Rt\ where j / k, as distinct if 

> 0.05. Approximate values of the CGLE parameters {R,v,^) at which any 
additional symmetry was gained or broken during the numerical continuation are listed in 
Table 3, only for those sequences where symmetry gaining or breaking behavior occurred. 

4.1 Features from the Solution Process 

To give an idea of the performance of the methodology employed, Figures 5-7 show several 
plots corresponding to application of the procedure for the sequence (see (11)), listed 
with id 8 in Tables 1-3. Continuation in this case was done on the CGLE parameter R 
only. Paths resulting from the continuation appear in Figure 5. Specifically, Figures 5a, 
5b, and 5c depict the resulting continuation paths by displaying, respectively, the values 
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id 

approximate {R,^,^) values: type of symmetry gained/broken 

2 

(16.5,-7,5): (16) gained —> (20.6,-7,5): (16) broken 
—> (53.4,-7, 5): (16) gained 

4 

(20.6,-7, 5): (16) broken —> (54.7,-7,5): (16) gained 

7 

(100,-7,0.3): (15), 1 = 2, gained —^ (100,-7,-0.02): (16), (17) gained 

8 

(32.6,-7, 5): (16) gained —> (85.5,-7,5): (16) broken 

10 

(82.8,-7,5.97): (16) broken —> (82.9,-7,5.97): (16) gained 
—^ (85.5,-7,5.97): (16) broken 

11 

(12.17,-7,5): (15), 1 = 2, (16), (17) gained —^ (60.1,-7, 5): (16), (17) broken 

13 

(16,-5.6, 3.4): (24) gained 

14 

(16.2,—7, 5): (16) gained —)• (15.8,—7,5): (15), 1 = 2, (17) gained 

—^ (80.6,-7, 5): (16), (17) broken —> (80.8,-7, 5): (16), (17) gained 
—^ (78.7,-7, 5): (16), (17) broken 

15 

(72.5,-7, 5): (16), (17) broken —^ (71.2,-7, 5): (16), (17) gained 

—^ (84.8,-7, 5): (16), (17) broken —> (67.7,-7, 5): (16), (17) gained 
—^ (80.9,-7, 5): (16), (17) broken —> (84.9,-7, 5): (16), (17) gained 


Table 3: Bifurcation points and type of symmetry gained/broken along the continuation 
paths. 

of the time period T, space translation S, and rotation® tp as functions of the continuation 
parameter R. Note from Figure 5b that within the range of i? ~ 32.6 through R ~ 85.5 
the value of S remained constant. The start of this interval of constant S corresponds to 
a step in the continuation process at which the resulting solution gained the additional 
symmetry (16); this symmetry was broken at the point in the path where S ceases to be 
constant. (Recall that solutions with symmetries (10) and (16) exist in subspaces of the 
solution space {A\ p, S, T) for which either 5 = 0 or 5 = T^/2; see Section 2.) 

The depictions in Figure 5 make it convenient to identify turning points in the continu¬ 
ation path and, for a given (fixed) value of the continuation parameter, whether there may 
exist multiple solutions of F = 0 in the path. For example, in Figure 5a one can identify 
four points where the line R = 90 intersects the curve T{R). These four points correspond 

®Since it was not strictly necessary, the constraint p G [0, 27r) was not explicitly enforced when solving 
the system of nonlinear algebraic equations (36). Furthermore, after performing a series of preliminary 
test runs, we found no advantage (from a computational point of view) in enforcing it. The values of tp in 
Figure 5c are displayed as they resulted from the solution of the system (36), and should be taken modulo 
an integer multiple of 27r, mapping them back to the interval [0, 27r). 
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0.098 



(a) T(R) 



Figure 5: Continuation sequence Paths traversed by T, S, and tp, as functions of R. 


to four invariant solutions computed at the same particular point (R, v, p) in parameter 
space. Then, the multitude of solutions associated with this point in parameter space can 
be inspected to determine whether they are different elements of the same orbit (41) of the 
symmetry group G, whether they belong to conjugate orbits of the symmetry group, or 
whether they belong to different (non-conjugate) orbits of the symmetry group. 

Spectra for several solutions in the path from i? = 16 to i? = 100 are shown in Figure 6. 
As expected, an increase in the value of R requires more terms in the expansions (12) and 
(34) in order to keep a suitable decay in both the spatial and temporal spectra for the 
solutions. Finally, surface plots of the real part 3^(A), imaginary part ^(A), and absolute 
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Figure 6: Spectra for three solutions in the sequence at values of i? = 16,60,100 
{u = -7, n = 5). 

value |j 4| for the solutions whose spectra are shown in Figure 6 appear in Figure 7, where 
the aforementioned gain and, thereafter, loss of symmetry (16) can be observed. 

The example above illustrates the general situation that one faces. By this we mean 
that, due to the use of the arc-length continuation option from the LOCA package [46], 
which was the appropriate choice for us because it allows for turning points in the path 
following process, it is possible (i.e., inherent in the continuation algorithm) that a point 
Pi = in the CGLE parameter space may return to itself, that is, p^+i = Pi, after 

k continuation steps. In such a situation, we must consider different cases for the solutions 
{\^^^■,ip{Pk+l)■,S{pk+l)lT{p,,+l)) and {Aj,^;ip{pi), S{pi),T{pi)). Namely, whether said solu¬ 
tions represent (i) the same G-orbit, that is, {<piPk+i), S{pk+i),T{pk+i)) = {(pipi), S{p,), T{pi)) 
and there exists some (0, a,T) £ G such that = {0, cr, t) ■ (ii) different, but conju¬ 
gate, G-orbits, as defined in Section 2 (see also (43)-(44)); or (iii) different, non-conjugate, 
G-orbits. Along a continuation path, say A^''\ many returns to a same point do occur. 
However, we include only one of the invariant solutions computed at pf' in the count 
in Table 2, since the presentation of the complete analysis of the multitude of invariant 
solutions that were computed at such “revisited” points is out of the scope of this paper. 

Challenging behavior that arose during the numerical continuation was often due to 
traversal of values of the continuation parameter in a cyclic manner, specifically related to 
the cases (i) and (ii) listed in the previous paragraph. As a result, the continuation path 
within these cycles would contain (different) elements in the same G-orbit, or solutions 
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(a) K(A), R=16 



X 

(b) 9(4), i? = 16 


X 

(c) |4|, R — 16 




(g) K(4), R = 100 



(h) 9(4), i?= 100 



(i) |4|, R= 100 


Figure 7; Surface plots of the real part 3^(A), imaginary part and absolute value 1^41 

for three solutions in sequence at values of v = —7, r = 5 and R = 16 (top), 72 = 60 
(middle), and R = 100 (bottom). Symmetry (16) was gained at {R,i',r) « (32.6,—7, 5) 
and broken at (72, i/,/i) (85.5,—7, 5). Hence A is even (about x = tt) for 72 = 60, but 

not for 72 = 16,100. Click here for a movie depicting the continuation path in (<y9, S, T) 
space, as well as solutions represented by plotting 9‘(H(x,0)) vs. )12(H(x,0)) at a sequence 
of continuation steps. 


representing conjugate G-orbits (cf. (14) and (43)-(44)). Within the cycles, solutions at 
the turning points in the continuation path were in the vicinity of solutions with additional 
symmetries, or near solutions with a smaller spatial period L^/q-i for some integer > 
1 or smaller time period Tjq^ for some integer ga > 1- Often the LOCA continuation 
algorithm [46] would exit from the cycles automatically, so that the procedure would again 
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start yielding solutions in different, non-conjugate, G-orbits, as well as continue to make 
progress towards the goal of reaching the (desired) final point in the CGLE parameter 
space. However, sometimes the continuation algorithm would get caught in said cycles. We 
discuss instances of these scenarios in the following paragraphs. 

An example of such cyclic behavior is depicted in Figure 8 for the sequence listed 

with id 11 in Tables 1-3, where continuation was again done on the CGLE parameter 
R. Figure 8a displays the values of the continuation parameter R, time period T, space 
translation S', and rotation ip as functions of the continuation step number. Traversal of 
repeated values for R, T, S, and ip is observed from the sub-figures in Figure 8a, where it 
can also be seen that the cycling behavior stops when i? ~ 12.17, at around continuation 
step number 180 (where S becomes constant), at which point the additional symmetries 
(15)-(17) are gained. 

Looking at Figure 8b, where the value of the space translation S is plotted as a func¬ 
tion of the continuation parameter R, one can see the cyclic behavior of R resulting in a 
symmetric curve with respect to the horizontal line at the vertical axis value of L^l2 = tt. 
The path depicted in Figure 8b, represented by the curve S{R), contains conjugate solu¬ 
tions (cf. (14) and (43)-(44)). More precisely, the numerical continuation path for values 
of R € [12, 20] that contains conjugate solutions is the one that yields the symmetric curve 
about the horizontal line at the value of L^l2 = tt (seen in Figure 8b). That is, points on 
the curve S(i?) that are mirror images with respect to the line Lj,/2 = vr correspond to 
conjugate solutions under the spatial reflection symmetry of the GGLE, which belong to 
conjugate orbits of the symmetry group G. The additional symmetry (16) was gained at a 
value of i? ~ 12.17, and at this point the cycling behavior stops and the spatial translation 
S takes on the value of L,,,/2, as solutions with symmetries (10) and (16) exist in subspaces 
of the solution space (A; y?, 5, T) for which either S' = 0 or 5 = L^/2 (refer to Section 2). 

Also, the additional symmetry (15), for I = 2, was gained along with the additional 
symmetry (16). Recall from (18) that the Fourier coefficients {am,n} of solutions with 
symmetry (15), for I = 2, satisfy am,n = 0 if m is even. Thus, we can visualize gain of 
this additional symmetry by selecting a coefficient am,n, for some even m and some n, and 
plotting its value as a function of the continuation parameter R, as done in Figure 8c for the 
coefficient Oq.o. At the point when this additional symmetry is gained, for R ~ 12.17, we see 
that the real part of the coefficient ao,o goes from (around) 0.22 to 0, whereas its imaginary 
part goes from (around) 0.1 to 0. Upon gaining the additional symmetry (15), for I = 2, 
the coefficient ao,o remains equal to zero, as seen in the path depicted in Figure 8c. Finally, 
from (18), a solution with symmetry (15), for I = 2, will have nonzero Fourier coefficients 
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(a) R, T, S, vs. continuation step number (b) S{R) 



(c) 5R(ao,o) vs. i?; ^(ao.o) vs. R 


(d) JR(a_i,o) vs. R-, ^(a-i^o) vs. R 


Figure 8: Cycling behavior during the numerical continuation for sequence (a) Values 

of the continuation parameter R and computed generator (ip, S, T) as functions of the 
continuation step number, (b) S' as a function of R. (c) 51?(ao,o) and ^(ao^o) as functions 
of R. (d) 51?(a_i,o) and 9(a_i^o) as functions of R. Click here for a movie depicting the 
continuation path in {cp, S, T) space, as well as solutions represented by plotting 0)) 

vs. 3^(A(x, 0)) at a sequence of continuation steps. The cycle in sub-figure (b) above can be 
seen during the initial steps of the path in {ip, S, T) space. Symmetries (15)-(17) were gained 
at {R,v,p) Ri (12.17,-7,5); symmetries (16)-(17) were broken at {R,v,r) k, (60.1,—7, 5). 


{am,n} for odd m. This is depicted for the coefficient a_i^ 0 ) plotted as a function of the 
continuation parameter R, in Figure 8d. As seen, the coefficient d_i,o remains nonzero after 
the additional symmetry (15) is gained (at the same time when the cycling behavior stops) 
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at a value of i? ~ 12.17. 

The aforementioned traversal of values of the continuation parameter in a cyclic manner 
was quite common, and often the LOCA continuation algorithm [46] would exit from the 
cycles automatically, that is, without us having to stop and restart the continuation with 
different values for the allowed increments on the continuation parameter. Nevertheless, as 
an alternative for circumventing such cycling behavior, we also experimented with taking 
the other parameters // or i/ in the CGLE (1) as continuation parameters. The continuation 
was always done on a single parameter at a time, while still all solutions were numerically 
continued from the regime with parameter values = (16,—7, 5) to the regime for 

= (100,-7,5). As a starting point for performing continuation on an alternate 
parameter, we would select a solution within the cycle for which the spectra (spatial or 
temporal, as appropriate) did not display characteristics typical of that of solutions around 
the turning points in the cycle. As an example, with the solutions represented via (34), 
given an integer > 1, the Fourier coefficients of a solution with spatial period 

L^/qi have a recognizable pattern of zeros, namely, = 0 if m is not divisible by qi. 
So if the numerical continuation was caught in a cycle where solutions around a turning 
point were close to a solution with spatial period L^/qi, as a starting point for performing 
continuation on an alternate parameter we could select a solution within the cycle for which 
Yhn > e, for m = 0, ±1,..., ±qi, and some cutoff, say, e = 10“^. 

One particular case in which it was beneficial to alternate the continuation parameter 
was for the sequence for which exiting automatically from cycling behavior in the 

vicinity of a solution that was even and had space period L^/2 and time period of T/2 
was challenging. Hence we experimented with alternating the continuation parameter, as 
indicated in the previous paragraph. The continuation then led to a solution with the 
additional symmetry (24), along with the invariance (10). This additional symmetry was 
gained at parameter values of i? = 16, ~ —5.6,/x ~ 3.4. Patterns resulting from the 

additional symmetry (24) can be visualized from the surface plot shown in Figure 9. 

4.2 Comments on Numerical Aspects 

In the present study, values of G [48,128] and Nt G [48,128] were used, respectively, in 
the truncation of the spatial Fourier series expansion (12) and the representation (34). (For 
comparison, values of = 32 and Nt = 48, 64 were used in the preceding study [35].) The 
number of terms Nx, Nt in each expansion was chosen so that the solutions had a decay of 
at least 10“® in their spatial and temporal spectra. The resulting number of unknowns for 
the system F = 0 of nonlinear algebraic equations (35) ranged between 4,000 and 32,260. 
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Figure 9: Surface plot of the absolute value |^| for a solution in sequence at parameter 
values = (100,-7,5). The plot is over the space period [O,^^,] on the horizontal 

axis and time period [0, T] on the vertical axis. The solution has the additional symmetry 
(24), with c = L^,(p = (p/2 + 7r,T = Tl2, so \A{x^ t)\ = \A{L^ — x,t + T/2)\, as observed 
from the plot. Click here for a movie depicting the continuation path in {(p, S, T) space, as 
well as surface plots of |^| at a sequence of continuation steps. The additional symmetry 
(24) was gained at k, (16,—5.6, 3.4). 



To solve the linear systems in (50) (related to the construction of the matrix Jj^Jr 
and vector in (49)) using the GMRES iterative solver from the Meschach library [47], 
we set a tolerance of 10“® for the residual || JsZfe — 6 fe ||2 and a maximum of 3,000 GMRES 
iterations. (Recall that the solution of said linear systems is needed for the computation 
of the Newton step, as defined by (49).) The number of iterations taken by the GMRES 
solver to meet the specihed residual tolerance ranged between 90 and 2,700. In terms of 
actual computing time (on a ThinkPad W530 personal workstation with 2.70 GHz processor 
speed), this translated to fractions of a second on the lower end to around 45-60 seconds on 
the higher end for the total time taken to compute the Newton step using (49). Occasionally 
the maximum number of GMRES iterations was reached for (a subset of) the four linear 
systems in (50) (related to the construction of the matrix Jj^Jr and vector J~^b in (49)). In 
this case the computation of the Newton step was reported as failed to the main numerical 
continuation routine. However, in most cases convergence to the desired residual tolerance 
was reached with under 2,000 GMRES iterations. Solution of the system of nonlinear 
algebraic equations typically took 2~6 iterations for Newton’s method (a maximum of 10 
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Newton iterations was set). Upon convergence of Newton’s method, the residual ||F ||2 was 
on the order of 10“^ or less. 

As for the LOCA numerical continuation library [46], we performed single-parameter 
continuation using the option of arc-length continuation in order to allow for turning points 
in the path following process. Other input information required by the LOCA library 
was set based on behavior observed for some initial runs. In particular, we experimented 
with providing the LOCA library values in the range [—0.1,1.0] for the initial change in 
the continuation parameter and [0.5,2.0] for the maximum increment in the continuation 
parameter. We found that it was best to set the initial change in the continuation parameter 
to be in the range [0.01,0.05], and to allow a maximum increment in the range [0.5,1.0]. 
Although larger values could also perform satisfactorily, in general we found that it was 
best for our problem to keep somewhat tight control on these increments in the sense that 
the number of failed attempts was then minimal (often zero). In addition, allowing large 
increments led several of the solutions to a single-frequency solution in the range R G [9,10] 
of values of the continuation parameter. With tighter bounds on the allowed increments, 
the numerical continuation led to a larger variety of solutions, as discussed at the beginning 
of Section 4 and in Section 4.1. 

Upon reaching a solution of F = 0 at the final CGLE parameter values, the values of 
Nx and Nt in the truncated expansions (12) and (34) were increased to confirm that, with 
the increased number of terms in the expansions, Newton’s method would converge to the 
same solution. (In other words, to confirm that the solution of F = 0 was numerically well 
defined.) In addition, the solution was validated against time integration of the truncated 
system of ODEs (28). Finally, we note that the computations were performed on a Thinkpad 
W530 personal workstation with 16 GB memory, four cores, with two threads per core, and 
2.70 GHz processor speed. Per the discussion in Section 3.2.1, four threads were used 
concurrently when solving for the Newton step using (49). 

5 Summary and Concluding Remarks 

We have considered the problem of computing invariant solutions of evolutionary partial 
differential equations (PDFs) with continuous symmetries via numerical continuation. We 
worked in particular with the complex Ginzburg-Landau equation (GGLE) 

BA B'^ A 

— = 72A + (l + m)^-(l + iM)A|A|2 
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in 1+1 space-time dimension with periodic boundary conditions in space, which has a three- 
parameter group of continuous symmetries generated by space-time translations x —?• x + a, 
t —>■ t + r and a rotation A —>■ A of the complex field A, in addition to being invariant 
under the action of the discrete group of transformations A{x^ t) —)• A{—x, t) of spatial 
reflections. We note, however, that it should be straightforward to apply the methodology 
employed to other evolutionary parameter-dependent PDEs invariant under the action of 
a group of continuous transformations. 

The desired solutions were sought from those of an underdetermined system of nonlinear 
algebraic equations resulting from the use of a spectral-Galerkin discretization of the CGLE 
in both space and time. Given a set of previously computed unstable invariant solutions 
of the GGLE for fixed values of the parameters R, v, and jjL [35], we employed numerical 
continuation to carry solutions of this initial set into solutions in a regime obtained by 
increasing the value of the parameter R. In other words, we considered the use of a path 
following, or homotopy, method to solve the system of nonlinear algebraic equations, with 
the initial guesses for the nonlinear equations solver (Newton’s method) provided by a set 
of solutions from [35]. We note that the GGLE exhibits chaotic behavior both at the initial 
and final parameter regions. 

The set of initial (unstable) solutions led to distinct, new invariant solutions of the 
CGLE in the final (and intermediate) parameter region(s). The resulting solutions are 
unstable, and have multiple modes and frequencies active in their spatial and temporal 
spectra, respectively. Thus this work and resulting solutions are an addition to the more 
extensively studied single-frequency and generalized traveling waves solutions of the CGLE. 
Symmetry gaining and breaking behavior, associated mainly with the spatial reflection 
symmetry A{x, t) —)> A{—x, t) of the CGLE, was detected frequently in the parameter 
regions traversed. This phenomenon is tantamount to occurrence of “phase transitions” 
in the space of G-orbits of solutions of the CGLE at particular values of the parameters 
(i?, +), and is one interesting aspect to consider as a focus for a subsequent detailed 

investigation. 

As discussed in Section 3.2.1, the use of conceptually simple, but effective techniques, 
led to an efficient computation of the Newton step. Employing POSIX threads (pthreads) 
programming [36] allowed us to take advantage of the multiple cores available nowadays 
in personal workstations and made the solution of a computationally challenging problem 
with a large number of unknowns (up to 32,260 unknowns in the present study) practical 
without the need of a cluster or supercomputer. Freely available software, specifically the 
Library of Continuation Algorithms (LOCA) package [46], the Meschach library [47] for 
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solving systems of linear algebraic equations, and the FFTW library [20] for computing 
fast Fourier transforms, were used as tools in the implementation of the software required 
to carry out the numerical continuation. 

Among aspects for further consideration we mention research on techniques that may 
help in minimizing or circumventing traversal of parameter values in a cyclic manner, per 
the discussion in Section 4.1. This could include alternative techniques for control of the 
step size in the continuation parameter or the use of multi-parameter continuation. Such 
additional features will provide a more versatile setting in which to explore further larger 
parameter regions (with an increasing number of unknowns and/or higher space dimension), 
and enhance the understanding of the structure of the solution space of the CGLE, and in 
particular, the structure of the space of orbits of its symmetry group. In addition, the fact 
that the resulting solutions are unstable suggests that the solutions may belong to the set 
of (infinitely many) unstable periodic orbits embedded in chaotic attractors [12, 10, 35]. 
This direction, by itself, is certainly very interesting to pursue in further studies of the 
dynamics of the CGLE, and on the potential use of such periodic orbits in the study of 
chaotic dynamical systems [12, 11, 55, 10]. 
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